##########################################
# Replication Data for Proksch, Lowe, Wäckerle, Soroka. (2018). Multilingual Sentiment Analysis: A New Approach to Measuring Conflict in Legislative Speeches. Legislative Studies Quarterly, Forthcoming.
##########################################

#Part 3: Wordscores Replication
#Most of this code follows the replication provided by Herzog and Benoit for the article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"
#We do not replicate the selection and Heckman models.

rm(list = ls(all = TRUE))
library(rstudioapi)

current_path <- getActiveDocumentContext()$path 
setwd(dirname(current_path ))
library(arm)
library(lme4)
library(sampleSelection)

### PATHS ########################################
dataAll <- "3_working_data_all.RData"
dataAll_senti <- "3_working_data_all_senti.RData"
dataSub <- "3_working_data_speakers.RData"
dataSub_senti <- "3_working_data_speakers_senti.RData"
dataOut <- "3_all_models_ML.RData"
##################################################

load(dataAll)
load(dataSub)
load(dataAll_senti)
load(dataSub_senti)


# Outcome models
# --------------
print("Estimating outcome models")    

pos1.ml <- lmer(data=dataSub, textscore ~ unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + (1 | m), control=lmerControl(optimizer = "bobyqa"))
summary(pos1.ml)
pos2.ml <- lmer(data=dataSub, textscore ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + government + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))
summary(pos2.ml)

pos3.ml <- lmer(data=dataSub, textscore ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))

# Outcome models with Sentiment
# --------------
print("Estimating outcome models with sentiment")    

pos1.senti.ml <- lmer(data=dataSub_senti, Sentiment ~ unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + (1 | m), control=lmerControl(optimizer = "bobyqa"))
summary(pos1.senti.ml)
pos2.senti.ml <- lmer(data=dataSub_senti, Sentiment ~ crisis*unemployment + crisis*safetyScaled + seniorityYearsScaled + government + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))

pos3.senti.ml <- lmer(data=dataSub_senti, Sentiment ~ crisis*government*unemployment + crisis*government*safetyScaled + seniorityYearsScaled + cabMember + (1 | m), control=lmerControl(optimizer = "bobyqa"))


save(pos1.ml,
     pos2.ml,
     pos3.ml,
     pos1.senti.ml,
     pos2.senti.ml,
     pos3.senti.ml,
     file=dataOut)

library(rjags)
set.seed(42)

### PATHS ########################################
dataSub <- "3_working_data_speakers.RData"
dataSub_senti <- "3_working_data_speakers_senti.RData"

##################################################

# MCMC values
mcmc.adapt <- 500
mcmc.burn <- 1000
mcmc.sample <- 2000
mcmc.thin <- 2
mcmc.chains <- 2

# JAGS model 
JagsModel <- "model{
for (i in 1:N)
{
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- a[m[i]] + inprod(b, X[i,1:K])
}

for (j in 1:J)
{
  a[j] ~ dnorm(mu.a, tau.a)
}

for (k in 1:K)
{
  b[k] ~ dnorm(0, 0.0001)
}

tau <- pow(sigma, -2)
sigma ~ dunif(0, 1000)

mu.a ~ dnorm(0, 0.0001)
tau.a <- pow(sigma.a, -2)
sigma.a ~ dunif(0, 1000)


}"


# function for JAGS data for multilevel models
forJagsData <- function(DV, xMatrix, groupID) {
  list(
    y = DV,
    X = xMatrix,
    N = nrow(xMatrix), 
    K = ncol(xMatrix),
    J = length(unique(groupID)),
    m = groupID
  )
}

load(dataSub)
data <- dataSub
load(dataSub_senti)
data_senti <- dataSub_senti


# Model 1 without interaction effects
# -----------------------------------
print("Estimating outcome model 1")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis
X <- model.matrix(eq, data)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data$textscore, X, data$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos1.posterior <- coda.samples(model,
                               variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                               n.iter = mcmc.sample,
                               thin = mcmc.thin)

save(pos1.posterior,file="3_outcome_model1_posterior.RData")


# Model 1 without interaction effects - SENTIMENT
# -----------------------------------
print("Estimating outcome model 1 with sentiment")

eq <- Sentiment ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis
X <- model.matrix(eq, data_senti)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data_senti$Sentiment, X, data_senti$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos1.senti.posterior <- coda.samples(model,
                                     variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                     n.iter = mcmc.sample,
                                     thin = mcmc.thin)

save(pos1.senti.posterior,file="3_outcome_model1_senti_posterior.RData")



# Model 2 with interaction effects
# ---------------------------------
print("Estimating outcome model 2")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data$textscore, X, data$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos2.posterior <- coda.samples(model,
                               variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                               n.iter = mcmc.sample,
                               thin = mcmc.thin)

save(pos2.posterior,file="3_outcome_model2_posterior.RData")

# Model 2 with interaction effects - SENTIMENT
# ---------------------------------
print("Estimating outcome model 2 with sentiment")

eq <- Sentiment ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled

X <- model.matrix(eq, data_senti)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data_senti$Sentiment, X, data_senti$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos2.senti.posterior <- coda.samples(model,
                                     variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                     n.iter = mcmc.sample,
                                     thin = mcmc.thin)

save(pos2.senti.posterior,file="3_outcome_model2_senti_posterior.RData")



# Model 3 with three-way interaction effects
# ------------------------------------------
print("Estimating outcome model 3")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled + crisis*government + government*unemployment + government*safetyScaled + crisis*government*unemployment + crisis*government*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data$textscore, X, data$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos3.posterior <- coda.samples(model,
                               variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                               n.iter = mcmc.sample,
                               thin = mcmc.thin)

save(pos3.posterior,file="3_outcome_model3_posterior.RData")

# Model 3 with three-way interaction effects - SENTIMENT
# ------------------------------------------
print("Estimating outcome model 3 with sentiment")

eq <- Sentiment ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled + crisis*government + government*unemployment + government*safetyScaled + crisis*government*unemployment + crisis*government*safetyScaled

X <- model.matrix(eq, data_senti)

model <- jags.model(
  textConnection(JagsModel),
  data = forJagsData(data_senti$Sentiment, X, data_senti$m),
  n.chains = mcmc.chains,
  n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos3.senti.posterior <- coda.samples(model,
                                     variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                     n.iter = mcmc.sample,
                                     thin = mcmc.thin)

save(pos3.senti.posterior,file="3_outcome_model3_senti_posterior.RData")


